What Was Tested: Do crop yields increase from GCOs globally?
What This Allows Me To Do: Identify if yield-distance relationships (near analysis) are a function of evolutionary escape or other variables.
What Are The Model Assumptions: GLS models are used as they have less assumptions about error structure and heteroscadicity.
#Barley -
#Select Data and Rescale Distance
near<-read.csv(file = "Near_All/Barley_Near_Centroid.csv", header=T)
near$rescale_ND <- near$NEAR_DIST/(1000*1000)
barley.near<-near %>% select(mean_barle,rescale_ND)
barley.near<-barley.near %>% filter(mean_barle > 0, na.rm=TRUE)
barley.near<-barley.near %>% filter(rescale_ND > 0, na.rm=TRUE)
barley.near$logBarley<-log10(barley.near$mean_barle)
#Plot
#Plot
ggplot(barley.near, aes(x=rescale_ND, y=logBarley)) +
geom_point(alpha=0.1, color="#00A9B4") + geom_smooth(method = "lm",color="#DB3500") + labs(color='GCO') +theme_justin +xlab("Distance From GCO (1000's km)") +ylab ("Log10 Yield")
## `geom_smooth()` using formula 'y ~ x'
Barley Near near Model,
### Model - Barley
# regular OLS no variance structure
barley.ols <- gls(logBarley ~ rescale_ND, data = barley.near)
# varFixed (variance changes linearly with X)
barley.fixed <- update(barley.ols, .~., weights = varFixed(~rescale_ND))
# varPower (variance changes as a power function with X)
barley.power <- update(barley.ols, . ~ ., weights = varPower(form = ~rescale_ND))
# varExp (variance changes as an exponential function of x)
barley.exp <- update(barley.ols, . ~., weights = varExp(form = ~ rescale_ND)) # error
# varConstPower (constant plus a power function of X (useful if X includes 0))
barley.ConstPower <- update(barley.ols, . ~., weights = varConstPower(form = ~ rescale_ND))
# compare all models by AIC
AIC(barley.ols, barley.fixed, barley.power, barley.ConstPower, barley.exp)
## df AIC
## barley.ols 3 19822.15
## barley.fixed 3 25277.66
## barley.power 4 19729.77
## barley.ConstPower 5 19731.77
## barley.exp 4 19730.85
#varPower
summary(barley.exp)
## Generalized least squares fit by REML
## Model: logBarley ~ rescale_ND
## Data: barley.near
## AIC BIC logLik
## 19730.85 19760.14 -9861.423
##
## Variance function:
## Structure: Exponential of variance covariate
## Formula: ~rescale_ND
## Parameter estimates:
## expon
## -0.01694584
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) -0.1452756 0.01100551 -13.20026 0
## rescale_ND 0.0147758 0.00138106 10.69889 0
##
## Correlation:
## (Intr)
## rescale_ND -0.867
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -4.5651801 -0.4496173 0.2959025 0.7252327 1.5929104
##
## Residual standard error: 0.6499998
## Degrees of freedom: 11197 total; 11195 residual
#Cassava -
#Select Data and Rescale Distance
near<-read.csv(file = "Near_All/Cassava_Near_Centroid.csv", header=T)
near$rescale_ND <- near$NEAR_DIST/(1000*1000)
cassava.near<-near %>% select(mean_cassa,rescale_ND)
cassava.near<-cassava.near %>% filter(mean_cassa > 0, na.rm=TRUE)
cassava.near<-cassava.near %>% filter(rescale_ND > 0, na.rm=TRUE)
cassava.near$logCassava<-log10(cassava.near$mean_cassa)
summary(cassava.near)
## mean_cassa rescale_ND logCassava
## Min. : 0.01222 Min. : 0.02601 Min. :-1.9128
## 1st Qu.: 3.67075 1st Qu.: 3.44972 1st Qu.: 0.5648
## Median : 7.73612 Median :10.39709 Median : 0.8885
## Mean : 8.17809 Mean :10.58149 Mean : 0.7396
## 3rd Qu.:12.24387 3rd Qu.:16.97533 3rd Qu.: 1.0879
## Max. :38.89642 Max. :19.72692 Max. : 1.5899
#Plot
ggplot(cassava.near, aes(x=rescale_ND, y=logCassava)) +
geom_point(alpha=0.2, color="#00A9B4") + geom_smooth(method = "lm",color="#DB3500") +theme_justin +xlab("Distance From GCO (1000's km)") +ylab ("Log10 Yield")
## `geom_smooth()` using formula 'y ~ x'
Cassava Near near Model,
### Model - Cassava
# regular OLS no variance structure
cassava.ols <- gls(logCassava ~ rescale_ND, data = cassava.near)
# varFixed (variance changes linearly with X)
cassava.fixed <- update(cassava.ols, .~., weights = varFixed(~rescale_ND))
# varPower (variance changes as a power function with X)
cassava.power <- update(cassava.ols, . ~ ., weights = varPower(form = ~rescale_ND))
# varExp (variance changes as an exponential function of x)
cassava.exp <- update(cassava.ols, . ~., weights = varExp(form = ~ rescale_ND)) # error
# varConstPower (constant plus a power function of X (useful if X includes 0))
cassava.ConstPower <- update(cassava.ols, . ~., weights = varConstPower(form = ~ rescale_ND))
# compare all models by AIC
AIC(cassava.ols, cassava.fixed, cassava.power, cassava.ConstPower, cassava.exp)
## df AIC
## cassava.ols 3 8900.148
## cassava.fixed 3 13094.247
## cassava.power 4 8862.351
## cassava.ConstPower 5 8864.351
## cassava.exp 4 8756.726
#varExp
summary(cassava.exp)
## Generalized least squares fit by REML
## Model: logCassava ~ rescale_ND
## Data: cassava.near
## AIC BIC logLik
## 8756.726 8783.727 -4374.363
##
## Variance function:
## Structure: Exponential of variance covariate
## Formula: ~rescale_ND
## Parameter estimates:
## expon
## -0.02047272
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 0.5155451 0.013402171 38.46728 0
## rescale_ND 0.0207906 0.000983958 21.12954 0
##
## Correlation:
## (Intr)
## rescale_ND -0.895
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -5.5684850 -0.3568539 0.3005624 0.7105753 1.7183862
##
## Residual standard error: 0.5998532
## Degrees of freedom: 6315 total; 6313 residual
#Groundnut -
#Select Data and Rescale Distance
near<-read.csv(file = "Near_All/Groundnut_Near_Centroid.csv", header=T)
near$rescale_ND <- near$NEAR_DIST/(1000*1000)
groundnut.near<-near %>% select(mean_groun,rescale_ND)
groundnut.near<-groundnut.near %>% filter(mean_groun > 0, na.rm=TRUE)
groundnut.near<-groundnut.near %>% filter(rescale_ND > 0, na.rm=TRUE)
groundnut.near$logGroundnut<-log10(groundnut.near$mean_groun)
#Plot
ggplot(groundnut.near, aes(x=rescale_ND, y=logGroundnut)) +
geom_point(alpha=0.2, color="#00A9B4") + geom_smooth(method = "lm",color="#DB3500") +theme_justin +xlab("Distance From GCO (1000's km)") +ylab ("Log10 Yield")
## `geom_smooth()` using formula 'y ~ x'
Groundnut Near near Model,
### Model - Groundnut
# regular OLS no variance structure
groundnut.ols <- gls(logGroundnut ~ rescale_ND, data = groundnut.near)
# varFixed (variance changes linearly with X)
groundnut.fixed <- update(groundnut.ols, .~., weights = varFixed(~rescale_ND))
# varPower (variance changes as a power function with X)
groundnut.power <- update(groundnut.ols, . ~ ., weights = varPower(form = ~rescale_ND))
# varExp (variance changes as an exponential function of x)
groundnut.exp <- update(groundnut.ols, . ~., weights = varExp(form = ~ rescale_ND)) # error
# varConstPower (constant plus a power function of X (useful if X includes 0))
groundnut.ConstPower <- update(groundnut.ols, . ~., weights = varConstPower(form = ~ rescale_ND))
# compare all models by AIC
AIC(groundnut.ols, groundnut.ConstPower,groundnut.power, groundnut.fixed, groundnut.exp)
## df AIC
## groundnut.ols 3 13520.01
## groundnut.ConstPower 5 13523.68
## groundnut.power 4 13521.68
## groundnut.fixed 3 18936.25
## groundnut.exp 4 13520.79
#varPower
summary(groundnut.exp)
## Generalized least squares fit by REML
## Model: logGroundnut ~ rescale_ND
## Data: groundnut.near
## AIC BIC logLik
## 13520.79 13548.9 -6756.397
##
## Variance function:
## Structure: Exponential of variance covariate
## Formula: ~rescale_ND
## Parameter estimates:
## expon
## 0.001516301
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) -0.24712201 0.013296847 -18.58501 0e+00
## rescale_ND 0.00429596 0.001069283 4.01761 1e-04
##
## Correlation:
## (Intr)
## rescale_ND -0.894
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -4.4934597 -0.4081239 0.2516889 0.6904358 1.8240400
##
## Residual standard error: 0.5352301
## Degrees of freedom: 8321 total; 8319 residual
#Maize -
#Select Data and Rescale Distance
near<-read.csv(file = "Near_All/Maize_Near_Centroid.csv", header=T)
near$rescale_ND <- near$NEAR_DIST/(1000*1000)
maize.near<-near %>% select(mean_maize,rescale_ND)
maize.near<-maize.near %>% filter(mean_maize > 0, na.rm=TRUE)
maize.near<-maize.near %>% filter(rescale_ND > 0, na.rm=TRUE)
maize.near$logMaize<-log10(maize.near$mean_maize)
#Plot
ggplot(maize.near, aes(x=rescale_ND, y=logMaize)) +geom_point(alpha=0.2, color="#F0AB00") +geom_smooth(method = "lm",color="#DB3500") +theme_justin +xlab("Distance From GCO (1000's km)") +ylab ("Log10 Yield")
## `geom_smooth()` using formula 'y ~ x'
Maize Near near Model,
### Model - Maize
# regular OLS no variance structure
maize.ols <- gls(logMaize ~ rescale_ND, data = maize.near)
# varFixed (variance changes linearly with X)
maize.fixed <- update(maize.ols, .~., weights = varFixed(~rescale_ND))
# varPower (variance changes as a power function with X)
maize.power <- update(maize.ols, . ~ ., weights = varPower(form = ~rescale_ND))
# varExp (variance changes as an exponential function of x) - Does not Converge
maize.exp <- update(maize.ols, . ~., weights = varExp(form = ~ rescale_ND)) # error
# varConstPower (constant plus a power function of X (useful if X includes 0))
maize.ConstPower <- update(maize.ols, . ~., weights = varConstPower(form = ~ rescale_ND))
# compare all models by AIC
AIC(maize.ols, maize.fixed, maize.power, maize.ConstPower,maize.exp)
## df AIC
## maize.ols 3 23829.13
## maize.fixed 3 26600.48
## maize.power 4 23786.72
## maize.ConstPower 5 23788.72
## maize.exp 4 23830.41
#varExp
summary(maize.exp)
## Generalized least squares fit by REML
## Model: logMaize ~ rescale_ND
## Data: maize.near
## AIC BIC logLik
## 23830.41 23860.54 -11911.2
##
## Variance function:
## Structure: Exponential of variance covariate
## Formula: ~rescale_ND
## Parameter estimates:
## expon
## 0.001294945
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 0.4023974 0.01228996 32.74195 0
## rescale_ND -0.0264575 0.00112441 -23.53007 0
##
## Correlation:
## (Intr)
## rescale_ND -0.918
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -4.6574695 -0.4741918 0.2375453 0.7011139 2.0342761
##
## Residual standard error: 0.5660735
## Degrees of freedom: 13792 total; 13790 residual
#Millet -
#Select Data and Rescale Distance
near<-read.csv(file="Near_All/Millet_Near_Centroid.csv", header=T)
near$rescale_ND <- near$NEAR_DIST/(1000*1000)
millet.near<-near %>% select(mean_mille,rescale_ND)
millet.near<-millet.near %>% filter(mean_mille > 0, na.rm=TRUE)
millet.near<-millet.near %>% filter(rescale_ND > 0, na.rm=TRUE)
millet.near$logMillet<-log10(millet.near$mean_mille)
#Plot
ggplot(millet.near, aes(x=rescale_ND, y=logMillet)) + geom_point(alpha=0.2, color="#FFB100") + geom_smooth(method = "lm",color="#DB3500") +theme_justin +xlab("Distance From GCO (1000's km)") +ylab ("Log10 Yield")
## `geom_smooth()` using formula 'y ~ x'
Millet Near near Model,
### Model - Millet
# regular OLS no variance structure
millet.ols <- gls(logMillet ~ rescale_ND, data = millet.near)
# varFixed (variance changes linearly with X)
millet.fixed <- update(millet.ols, .~., weights = varFixed(~rescale_ND))
# varPower (variance changes as a power function with X)
millet.power <- update(millet.ols, . ~ ., weights = varPower(form = ~rescale_ND))
# varExp (variance changes as an exponential function of x)
millet.exp <- update(millet.ols, . ~., weights = varExp(form = ~ rescale_ND)) # error
# varConstPower (constant plus a power function of X (useful if X includes 0))
millet.ConstPower <- update(millet.ols, . ~., weights = varConstPower(form = ~ rescale_ND))
# compare all models by AIC
AIC(millet.ols, millet.fixed, millet.power, millet.ConstPower, millet.exp)
## df AIC
## millet.ols 3 12847.29
## millet.fixed 3 15549.03
## millet.power 4 12843.80
## millet.ConstPower 5 12845.80
## millet.exp 4 12849.22
#varPower
summary(millet.power)
## Generalized least squares fit by REML
## Model: logMillet ~ rescale_ND
## Data: millet.near
## AIC BIC logLik
## 12843.8 12871.38 -6417.901
##
## Variance function:
## Structure: Power of variance covariate
## Formula: ~rescale_ND
## Parameter estimates:
## power
## 0.02436811
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) -0.3207978 0.011790772 -27.207533 0e+00
## rescale_ND -0.0060439 0.001533657 -3.940818 1e-04
##
## Correlation:
## (Intr)
## rescale_ND -0.816
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -4.9859199 -0.4021284 0.2580831 0.6758465 1.7350781
##
## Residual standard error: 0.560494
## Degrees of freedom: 7298 total; 7296 residual
#Rapeseed -
#Select Data and Rescale Distance
near<-read.csv(file="./Near_All/Rapeseed_Near_Centroid.csv", header=T)
near$rescale_ND <- near$NEAR_DIST/(1000*1000)
rapeseed.near<-near %>% select(mean_rapes,rescale_ND)
rapeseed.near<-rapeseed.near %>% filter(mean_rapes > 0, na.rm=TRUE)
rapeseed.near<-rapeseed.near %>% filter(rescale_ND > 0, na.rm=TRUE)
rapeseed.near$logRapeseed<-log10(rapeseed.near$mean_rapes)
#Plot
ggplot(rapeseed.near, aes(x=rescale_ND, y=logRapeseed)) +geom_point(alpha=0.2, color="#FFAF02") + geom_smooth(method = "lm",color="#DB3500") +theme_justin +xlab("Distance From GCO (1000's km)") +ylab ("Log10 Yield")
## `geom_smooth()` using formula 'y ~ x'
Rapeseed Near near Model,
### Model - Rapeseed
# regular OLS no variance structure
rapeseed.ols <- gls(logRapeseed ~ rescale_ND, data = rapeseed.near)
# varFixed (variance changes linearly with X)
rapeseed.fixed <- update(rapeseed.ols, .~., weights = varFixed(~rescale_ND))
# varPower (variance changes as a power function with X)
rapeseed.power <- update(rapeseed.ols, . ~ ., weights = varPower(form = ~rescale_ND))
# varExp (variance changes as an exponential function of x)
rapeseed.exp <- update(rapeseed.ols, . ~., weights = varExp(form = ~ rescale_ND)) # error
# varConstPower (constant plus a power function of X (useful if X includes 0))
rapeseed.ConstPower <- update(rapeseed.ols, . ~., weights = varConstPower(form = ~ rescale_ND))
# compare all models by AIC
AIC(rapeseed.ols, rapeseed.fixed, rapeseed.power, rapeseed.ConstPower, rapeseed.exp)
## df AIC
## rapeseed.ols 3 15008.97
## rapeseed.fixed 3 19957.49
## rapeseed.power 4 15008.34
## rapeseed.ConstPower 5 15010.34
## rapeseed.exp 4 15003.10
#varPower
summary(rapeseed.exp)
## Generalized least squares fit by REML
## Model: logRapeseed ~ rescale_ND
## Data: rapeseed.near
## AIC BIC logLik
## 15003.1 15031.64 -7497.552
##
## Variance function:
## Structure: Exponential of variance covariate
## Formula: ~rescale_ND
## Parameter estimates:
## expon
## -0.005850967
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) -0.24958856 0.010965742 -22.760754 0e+00
## rescale_ND -0.00563787 0.001513497 -3.725062 2e-04
##
## Correlation:
## (Intr)
## rescale_ND -0.857
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -4.5076091 -0.4305570 0.3144485 0.7080704 1.5508114
##
## Residual standard error: 0.5629072
## Degrees of freedom: 9258 total; 9256 residual
#Rice -
#Select Data and Rescale Distance
near<-read.csv(file="./Near_All/Rice_Near_Centroid.csv", header=T)
near$rescale_ND <- near$NEAR_DIST/(1000*1000)
rice.near<-near %>% select(mean_rice,rescale_ND)
rice.near<-rice.near %>% filter(mean_rice > 0, na.rm=TRUE)
rice.near<-rice.near %>% filter(rescale_ND > 0, na.rm=TRUE)
rice.near$logRice<-log10(rice.near$mean_rice)
#Plot
ggplot(rice.near, aes(x=rescale_ND, y=logRice)) +geom_point(alpha=0.2, color="#FBA900") + geom_smooth(method = "lm",color="#DB3500") +theme_justin +xlab("Distance From GCO (1000's km)") +ylab ("Log10 Yield")
## `geom_smooth()` using formula 'y ~ x'
Rice Near near Model,
### Model - Rice
# regular OLS no variance structure
rice.ols <- gls(logRice ~ rescale_ND, data = rice.near)
# varFixed (variance changes linearly with X)
rice.fixed <- update(rice.ols, .~., weights = varFixed(~rescale_ND))
# varPower (variance changes as a power function with X)
rice.power <- update(rice.ols, . ~ ., weights = varPower(form = ~rescale_ND))
# varExp (variance changes as an exponential function of x)
rice.exp <- update(rice.ols, . ~., weights = varExp(form = ~ rescale_ND)) # error
# varConstPower (constant plus a power function of X (useful if X includes 0))
rice.ConstPower <- update(rice.ols, . ~., weights = varConstPower(form = ~ rescale_ND))
# compare all models by AIC
AIC(rice.ols, rice.fixed, rice.power, rice.ConstPower, rice.exp)
## df AIC
## rice.ols 3 13556.55
## rice.fixed 3 14019.13
## rice.power 4 13242.15
## rice.ConstPower 5 13244.15
## rice.exp 4 13449.33
#varPower
summary(rice.power)
## Generalized least squares fit by REML
## Model: logRice ~ rescale_ND
## Data: rice.near
## AIC BIC logLik
## 13242.15 13270.31 -6617.077
##
## Variance function:
## Structure: Power of variance covariate
## Formula: ~rescale_ND
## Parameter estimates:
## power
## 0.1988213
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 0.3714575 0.008714763 42.62394 0
## rescale_ND -0.0291549 0.001030728 -28.28577 0
##
## Correlation:
## (Intr)
## rescale_ND -0.768
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -4.5128640 -0.4174723 0.1789334 0.6418429 1.7784249
##
## Residual standard error: 0.3694383
## Degrees of freedom: 8428 total; 8426 residual
#Rye -
#Select Data and Rescale Distance
near<-read.csv(file="./Near_All/Rye_Near_Centroid.csv", header=T)
near$rescale_ND <- near$NEAR_DIST/(1000*1000)
rye.near<-near %>% select(mean_rye_Y,rescale_ND)
rye.near<-rye.near %>% filter(mean_rye_Y > 0, na.rm=TRUE)
rye.near<-rye.near %>% filter(rescale_ND > 0, na.rm=TRUE)
rye.near$logRye<-log10(rye.near$mean_rye_Y)
#Plot
ggplot(rye.near, aes(x=rescale_ND, y=logRye)) +geom_point(alpha=0.2, color="#FFB200") +geom_smooth(method = "lm",color="#DB3500") +theme_justin +xlab("Distance From GCO (1000's km)") +ylab ("Log10 Yield")
## `geom_smooth()` using formula 'y ~ x'
Rye Near near Model,
### Model - Rye
# regular OLS no variance structure
rye.ols <- gls(logRye ~ rescale_ND, data = rye.near)
# varFixed (variance changes linearly with X)
rye.fixed <- update(rye.ols, .~., weights = varFixed(~rescale_ND))
# varPower (variance changes as a power function with X)
rye.power <- update(rye.ols, . ~ ., weights = varPower(form = ~rescale_ND))
# varExp (variance changes as an exponential function of x)
rye.exp <- update(rye.ols, . ~., weights = varExp(form = ~ rescale_ND)) # error
# varConstPower (constant plus a power function of X (useful if X includes 0))
rye.ConstPower <- update(rye.ols, . ~., weights = varConstPower(form = ~ rescale_ND))
# compare all models by AIC
AIC(rye.ols, rye.fixed, rye.power, rye.ConstPower, rye.exp)
## df AIC
## rye.ols 3 14495.61
## rye.fixed 3 18866.35
## rye.power 4 14483.81
## rye.ConstPower 5 14367.69
## rye.exp 4 14497.54
#varConstPower
summary(rye.ConstPower)
## Generalized least squares fit by REML
## Model: logRye ~ rescale_ND
## Data: rye.near
## AIC BIC logLik
## 14367.69 14402.87 -7178.846
##
## Variance function:
## Structure: Constant plus power of variance covariate
## Formula: ~rescale_ND
## Parameter estimates:
## const power
## 5.970222 -2.338874
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 0.02234046 0.010606468 2.106306 0.0352
## rescale_ND -0.03242745 0.001418423 -22.861619 0.0000
##
## Correlation:
## (Intr)
## rescale_ND -0.814
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -4.8648858 -0.5176796 0.2805052 0.6733118 2.1347298
##
## Residual standard error: 0.0924272
## Degrees of freedom: 8400 total; 8398 residual
#Sorghum -
#Select Data and Rescale Distance
near<-read.csv(file="./Near_All/Sorghum_Near_Centroid.csv", header=T)
near$rescale_ND <- near$NEAR_DIST/(1000*1000)
sorghum.near<-near %>% select(mean_sorgh,rescale_ND)
sorghum.near<-sorghum.near %>% filter(rescale_ND > 0, na.rm=TRUE)
sorghum.near<-sorghum.near %>% filter(mean_sorgh > 0, na.rm=TRUE)
sorghum.near$logSorghum<-log10(sorghum.near$mean_sorgh)
#Plot
ggplot(sorghum.near, aes(x=rescale_ND, y=logSorghum)) + geom_point(alpha=0.2, color="#00A9B4") +geom_smooth(method = "lm",color="#DB3500") +theme_justin +xlab("Distance From GCO (1000's km)") +ylab ("Log10 Yield")
## `geom_smooth()` using formula 'y ~ x'
Sorghum Near near Model,
### Model - Sorghum
# regular OLS no variance structure
sorghum.ols <- gls(logSorghum ~ rescale_ND, data = sorghum.near)
# varFixed (variance changes linearly with X)
sorghum.fixed <- update(sorghum.ols, .~., weights = varFixed(~rescale_ND))
# varPower (variance changes as a power function with X)
sorghum.power <- update(sorghum.ols, . ~ ., weights = varPower(form = ~rescale_ND))
# varExp (variance changes as an exponential function of x)
sorghum.exp <- update(sorghum.ols, . ~., weights = varExp(form = ~ rescale_ND)) # error
# varConstPower (constant plus a power function of X (useful if X includes 0))
sorghum.ConstPower <- update(sorghum.ols, . ~., weights = varConstPower(form = ~ rescale_ND))
# compare all models by AIC
AIC(sorghum.ols, sorghum.fixed, sorghum.power, sorghum.ConstPower, sorghum.exp)
## df AIC
## sorghum.ols 3 18017.89
## sorghum.fixed 3 22316.61
## sorghum.power 4 17964.38
## sorghum.ConstPower 5 17966.38
## sorghum.exp 4 17921.13
#Exp
summary(sorghum.exp)
## Generalized least squares fit by REML
## Model: logSorghum ~ rescale_ND
## Data: sorghum.near
## AIC BIC logLik
## 17921.13 17950.06 -8956.567
##
## Variance function:
## Structure: Exponential of variance covariate
## Formula: ~rescale_ND
## Parameter estimates:
## expon
## -0.01997319
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) -0.5346044 0.012882963 -41.49701 0
## rescale_ND 0.0535188 0.001519261 35.22687 0
##
## Correlation:
## (Intr)
## rescale_ND -0.897
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -4.9118332 -0.4762013 0.2568962 0.6495763 2.1113273
##
## Residual standard error: 0.6680175
## Degrees of freedom: 10229 total; 10227 residual
#Soybean -
#Select Data and Rescale Distance
near<-read.csv(file="./Near_All/Soybean_Near_Centroid.csv", header=T)
near$rescale_ND <- near$NEAR_DIST/(1000*1000)
soybean.near<-near %>% select(mean_soybe,rescale_ND)
soybean.near<-soybean.near %>% filter(mean_soybe > 0, na.rm=TRUE)
soybean.near<-soybean.near %>% filter(rescale_ND > 0, na.rm=TRUE)
soybean.near$logSoybean<-log10(soybean.near$mean_soybe)
#Plot
ggplot(soybean.near, aes(x=rescale_ND, y=logSoybean)) +geom_point(alpha=0.2, color="#00A9B4") +geom_smooth(method = "lm",color="#DB3500") +theme_justin +xlab("Distance From GCO (1000's km)") +ylab ("Log10 Yield")
## `geom_smooth()` using formula 'y ~ x'
Soybean Near near Model,
### Model - Soybean
# regular OLS no variance structure
soybean.ols <- gls(logSoybean ~ rescale_ND, data = soybean.near)
# varFixed (variance changes linearly with X)
soybean.fixed <- update(soybean.ols, .~., weights = varFixed(~rescale_ND))
# varPower (variance changes as a power function with X)
soybean.power <- update(soybean.ols, . ~ ., weights = varPower(form = ~rescale_ND))
# varExp (variance changes as an exponential function of x)
soybean.exp <- update(soybean.ols, . ~., weights = varExp(form = ~ rescale_ND)) # error
# varConstPower (constant plus a power function of X (useful if X includes 0))
soybean.ConstPower <- update(soybean.ols, . ~., weights = varConstPower(form = ~ rescale_ND))
# compare all models by AIC
AIC(soybean.ols, soybean.fixed, soybean.power, soybean.exp, soybean.ConstPower)
## df AIC
## soybean.ols 3 17503.72
## soybean.fixed 3 21913.61
## soybean.power 4 17499.79
## soybean.exp 4 17485.46
## soybean.ConstPower 5 17501.79
#varPower
summary(soybean.exp)
## Generalized least squares fit by REML
## Model: logSoybean ~ rescale_ND
## Data: soybean.near
## AIC BIC logLik
## 17485.46 17514.52 -8738.731
##
## Variance function:
## Structure: Exponential of variance covariate
## Formula: ~rescale_ND
## Parameter estimates:
## expon
## -0.006650274
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) -0.3391456 0.009434459 -35.94754 0
## rescale_ND 0.0164092 0.001066931 15.37981 0
##
## Correlation:
## (Intr)
## rescale_ND -0.821
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -5.4419479 -0.4105329 0.3369870 0.6747373 1.6719120
##
## Residual standard error: 0.5796952
## Degrees of freedom: 10548 total; 10546 residual
#Sunflower -
#Select Data and Rescale Distance
near<-read.csv(file="Near_All/Sunflower_Near_Centroid.csv", header=T)
near$rescale_ND <- near$NEAR_DIST/(1000*1000)
sunflower.near<-near %>% select(mean_sunfl,rescale_ND)
sunflower.near<-sunflower.near %>% filter(mean_sunfl > 0, na.rm=TRUE)
sunflower.near<-sunflower.near %>% filter(rescale_ND > 0, na.rm=TRUE)
sunflower.near$logSunflower<-log10(sunflower.near$mean_sunfl)
#Plot
ggplot(sunflower.near, aes(x=rescale_ND, y=logSunflower)) +geom_point(alpha=0.2, color="#FFB200") +geom_smooth(method = "lm",color="#DB3500") +theme_justin +xlab("Distance From GCO (1000's km)") +ylab ("Log10 Yield")
## `geom_smooth()` using formula 'y ~ x'
Sunflower Near near Model,
### Model - Sunflower
# regular OLS no variance structure
sunflower.ols <- gls(logSunflower ~ rescale_ND, data = sunflower.near)
# varFixed (variance changes linearly with X)
sunflower.fixed <- update(sunflower.ols, .~., weights = varFixed(~rescale_ND))
# varPower (variance changes as a power function with X)
sunflower.power <- update(sunflower.ols, . ~ ., weights = varPower(form = ~rescale_ND))
# varExp (variance changes as an exponential function of x)
sunflower.exp <- update(sunflower.ols, . ~., weights = varExp(form = ~ rescale_ND)) # error
# varConstPower (constant plus a power function of X (useful if X includes 0))
sunflower.ConstPower <- update(sunflower.ols, . ~., weights = varConstPower(form = ~ rescale_ND))
# compare all models by AIC
AIC(sunflower.ols, sunflower.fixed, sunflower.power,sunflower.ConstPower, sunflower.exp)
## df AIC
## sunflower.ols 3 15907.51
## sunflower.fixed 3 19567.89
## sunflower.power 4 15829.34
## sunflower.ConstPower 5 15831.34
## sunflower.exp 4 15881.72
#varPower
summary(sunflower.ConstPower)
## Generalized least squares fit by REML
## Model: logSunflower ~ rescale_ND
## Data: sunflower.near
## AIC BIC logLik
## 15831.34 15867.03 -7910.67
##
## Variance function:
## Structure: Constant plus power of variance covariate
## Formula: ~rescale_ND
## Parameter estimates:
## const power
## 5.417105e-06 8.663531e-02
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) -0.29574009 0.012994679 -22.758552 0
## rescale_ND -0.00778374 0.001392125 -5.591268 0
##
## Correlation:
## (Intr)
## rescale_ND -0.894
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -4.5510755 -0.4076590 0.2530632 0.7213569 1.6555183
##
## Residual standard error: 0.4766152
## Degrees of freedom: 9302 total; 9300 residual
#Wheat -
#Select Data and Rescale Distance
near<-read.csv(file="Near_All/Wheat_Near_Centroid.csv", header=T)
near$rescale_ND <- near$NEAR_DIST/(1000*1000)
wheat.near<-near %>% select(mean_wheat,rescale_ND)
wheat.near<-wheat.near %>% filter(mean_wheat > 0, na.rm=TRUE)
wheat.near<-wheat.near %>% filter(rescale_ND> 0, na.rm=TRUE)
wheat.near$logWheat<-log10(wheat.near$mean_wheat)
#Plot
ggplot(wheat.near, aes(x=rescale_ND, y=logWheat)) +
geom_point(alpha=0.2, color="#00A9B4") + geom_smooth(method = "lm",color="#DB3500")+theme_justin +xlab("Distance From GCO (1000's km)") +ylab ("Log10 Yield")
## `geom_smooth()` using formula 'y ~ x'
Wheat Near near Model,
### Model - Wheat
# regular OLS no variance structure
wheat.ols <- gls(logWheat ~ rescale_ND, data = wheat.near)
# varFixed (variance changes linearly with X)
wheat.fixed <- update(wheat.ols, .~., weights = varFixed(~rescale_ND))
# varPower (variance changes as a power function with X)
wheat.power <- update(wheat.ols, . ~ ., weights = varPower(form = ~rescale_ND))
# varExp (variance changes as an exponential function of x)
wheat.exp <- update(wheat.ols, . ~., weights = varExp(form = ~ rescale_ND)) # error
# varConstPower (constant plus a power function of X (useful if X includes 0))
wheat.ConstPower <- update(wheat.ols, . ~., weights = varConstPower(form = ~ rescale_ND))
# compare all models by AIC
AIC(wheat.ols, wheat.fixed, wheat.power, wheat.ConstPower, wheat.exp)
## df AIC
## wheat.ols 3 22293.30
## wheat.fixed 3 27927.21
## wheat.power 4 22226.63
## wheat.ConstPower 5 22228.63
## wheat.exp 4 22226.87
#varPower
summary(wheat.power)
## Generalized least squares fit by REML
## Model: logWheat ~ rescale_ND
## Data: wheat.near
## AIC BIC logLik
## 22226.63 22256.46 -11109.32
##
## Variance function:
## Structure: Power of variance covariate
## Formula: ~rescale_ND
## Parameter estimates:
## power
## -0.07538861
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) -0.10244907 0.010269171 -9.976371 0
## rescale_ND 0.00872257 0.001297664 6.721749 0
##
## Correlation:
## (Intr)
## rescale_ND -0.87
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -4.7594172 -0.4459564 0.2490299 0.6773372 1.7349711
##
## Residual standard error: 0.6526273
## Degrees of freedom: 12812 total; 12810 residual